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Abstract 

This paper surveys the state-of-the-art of turbulent flow 
computation. The formulations have been generalized somewhat to 
increase the range of their applicability, and the excitement of 
current debate on equation models has been brought into the 
review. Some new ideas on the modeling of the pressure-strain 
term in the Reynolds stress equations are also suggested. 

The review was prepared from lecture notes drafted for an 
A.I.Ch.E. short course on turbulence structure (Chicago, Nov., 
1970 ). 
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NOMENCLATURE 


% 
q / 2 


wall layer thickness parameter (4.9) 
structure tensor ( 6 . 11 ) 
isotropic dissipation ( 6 . 3 d) 
pressure gradient parameter ( 7 * 3 ) 
turbulence length scale 
turbulence production -R. .S. . 
pressure-strain tensor ( 8 . 1 ) 
mean pressure 
fluctuation pressure 
turbulence Prandtl number (4.20) 
pressure gradient parameter (4.10a) 
large eddy vector velocity scale 
turbulence kinetic energy density 
Reynolds number 

u.u. , "Reynolds stress tensor" 
mean strain rate tensor (4.2b) 
time 

fluctuation velocity vector (u,v,w) 

mean velocity vector (U,V, W) 

friction velocity* \/ p* 

transpiration parameter (4.10b) 

viscous terms in R. . equation ( 6 . 3a), ( 6 . 5a) 

-*-C . . 

cartesian coordinate vector (x,y, z) 

molectular thermal diffusivity 

turbulent thermal diffusivity 

molecular kinematic viscosity 

turbulent kinematic viscosity 

Karman constant 

mass density 

shear stress 

wall shear stress 

mean temperature 

fluctuation temperature 

see ( 6 . 5 b) 

dissipation of turbulence energy (6.4) 
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ACRONYMS 


CHC 

Champagne, Harris, and Corrsin 1970 


DH 

Daly and Harlow 1970 


DR 

Donaldson and Rosenbaum 1968 


FVF 

Fluctuating Velocity Field closure 


MH 

Mellor and Herring 1970a, b 


MVF 

Mean Velocity Field closure 


MVFN 

Newtonian MVF closure 


MTE 

Mean Turbulent Energy closure 


MTEN 

Newtonian MTE closure 


MTES 

Structural MTE closure 


MTEN/L 

MTEN closure with dynamical length scale 

equation 

MRS 

Mean Reynolds Stress closure 


MRS/L 

MRS closure with dynamical length scale 

equation 

TBLPC 

TBLPC I 968 


TR 

Tucker and Reynolds 1968 
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FIGURE CAPTIONS 


Fig, 1 
Fig. 2 
Pig. 3 

Fig. 4 
Pig. 5 
Fig. 6 
Fig. 7 

Fig. 8 

Fig. 9 

Fig. 10 

Fig. 11 

Fig. 12 

Fig. 13 

Fig. 14 

Fig. 15 
Fig. 16 

Fig. 17 

Fig. 18 


TBLPC test flow no. 1400 
TBLPC test flow no. 2600 

TBLPC evaluation committee rankings of the methods. 

A lower ranking is better. 

Wall layer thickness parameter as used by ^ays 
Typical mixing length distribution in boundary layers 
Patankar-Spalding MVFN wall jet prediction 
Dvorak MVFN calculation for a wall jet in a boundary 
layer subjected to strong adverse pressure gradient. 
Comparison of calculations on the symmetry plane in a 
three-dimensional flow (Wheeler and Johnston, unpublished) 
Herring and Mellor MVFN calculation for the skin friction 
factor on a flat plate in compressible flow 
Herring and Mellor MVFN calculation for a compressible 
boundary layer on a waisted body of revolution. a)integral 
parameters b) profiles 

Healzer and Kays (unpublished) MFVN calculation of the heat 
transfer coefficient in a supersonic rocket nozzle 
Seban and Chin MVFN calculation of the recirculating 
flow in a square cavity 

Hussain and Reynolds MVFN calculation of the dispersion 
relationship for plane waves in turbulent channel flow 
NCAR 6-layer atmospheric circulation model. a) model 
b) calculated sea level isobars 

Bradshaw's MTES calculation for the Tillmann Ledge flow 
Nash's MTES calculation of a three-dimensional turbulent 
boundary layer 

Mellor and Herring's MTEN calculation for a flat plate 
boundary layer a) inner region b) outer region 
Mellor and Herring's MTEN calculation for a relaxing 
boundary layer 
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Fig. 19 Kays' MTEN calculation for the heat transfer to an 
accelerating boundary layer 

Fig. 20 Kays 1 MTEN calculation for the heat transfer to an 
accelerating boundary layer with transpiration 
Fig. 21 Kays' MTEN calculation for the heat transfer to an 

accelerating boundary layer with changes in transpiration 
Fig. 22 Kearney's MTEN calculation of the effects of free stream 
turbulence on an accelerating boundary layer 
Fig. 23 Spalding and Rodi's MTEN calculation, incorporating their 
dynamical equation for the length scale, for a plane jet 
Fig. 24 Donaldson and Rosenbaum's MRS calculation for a flat 
plate boundary layer 

Fig. 25a Turbulence energy in the GHC flow-determination of C 2 
Fig. 25b Length scales in the CHC flow. Note that C 2 =2 models 
the observed length scale percentage changes 
Fig. 25c R^ . in the CHC flow; all use (9.1) and (9*6). A-(9-8) 
with C 3 =5. B- (9.9c) with C^=l/2. C- (9-9d) with 
C 4 =l/2, C 5 =l/4 

Fig. 25d Dissipation length scales in the CHC flow. A,B, C as for 
Fig. 25c. Note that C reproduces the behavior deduced 
using (9«1) and (9.6) 

Fig. 26a Turbulence energy in the TR flow-determination of C2 
Fig. 26b Structure in the unstrained return-to-isotropy portion 
of the TR flow; determination of C^ 

Fig. 26c Turbulence energy in the TR flow. A- (9.8) with C 3 =5. 

B- (9.9c) with C^=l/4 

Fig. 26d Structure in the straining region of the TR flow; A,B 
as for Fig. 26c. Note that B reproduces the length 
scale changes calculated using (9*1) and (9*6) 
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1 . Scope and Purpose of This Survey 


The objective of this survey is to provide a brief but 
reasonably complete account of the state-of-the-art of turbulent 
flow computations. The review will be limited to methods that 
have some scientific basis* that show promise for extension to 
wider classes of flows, and that have been developed to the point 
where at least some technical information of practical use can be 
obtained. Emphasis will be placed on the physical assumptions 
rather than on the numerical techniques. The central ideas of 
contemporary methods will be highlighted, and we shall direct the 
reader to individual sources for more detailed descriptions. An 
effort to both relate and critique the methods has been made. A 
significant portion of the material covered here is not yet pub- 
lished elsewhere, and this material should be valuable to workers 
in this field. 

Generators of computational schemes for turbulent flows are 
often most exuberant when scant data are available for comparison 
with their predictions. This is indeed the case for most classes 
of flows, with the single exception of two-dimensional steady in- 
compressible turbulent boundary layers. In I968 a turbulent 
boundary layer prediction method calibration conference (TBLPC 
I968) was held at Stanford, where for the first time a large 
number of methods (29) were compared on a systematic basis. This 
comparison established the viability of prediction methods based 
on various closure models of the partial differential equations 
describing turbulent boundary layer flows. Systematic extensions 
of such comparisons are now in progress in several quarters, and 
I am indebted to several colleagues for making their unpublished 
work available for this review. 

In other classes of flows, such as free shear layers, 
unsteady boundary layer flows, flows with strong boundary layer- 
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inviscid region interaction, with rotation, or buoyancy, separated 
flows, cavity flows, etc., the data available are very spotty, and 
our ability to evaluate such computations is therefore limited. 

But it does seem clear that there remains considerable room for 
improvement and extension of existing methods for these flows. 

Finally, there is increased activity in exploration of 
complex differential equation models. This effort has not been 
very systematic. Of particular concern is the so-called pressure- 
strain term in the dynamical equations for the Reynolds stresses. 

A systematic development of what may be a better model for this 
term is developed here for consideration by my colleagues. 
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2 . The Stanford Conference* 


The work leading up to the I 968 TBLPC produced a volume of 
target boundary layer data. A committee headed by D. Coles 
surveyed over 100 experiments and selected 33 flows for inclusion 
in this volume. The data of each experiment were carefully 
reanalyzed, recomputed for placement in a standard form, critiques 
were solicited from the experimentors, and all this was documented 
in a tidy manner by E. Hirst and D. Coles (TBLPC-2). These data 
now stand as a classic base of comparison for turbulent boundary 
layer prediction methods. Only the hydrodynamic aspects of these 
layers were considered, and a corresponding standard for thermal 
behavior is still lacking (though the data of W. M. Kays and his 
associates are rapidly becoming such a collection). ' Moreover, the 
flows selected were relatively mild. Very strong pressure gradi- 
ents, transpiration, roughness, rotations, and other interesting 
effects were not included. 

Sixteen flows were selected as mandatory computations (most 
predictors did the others as well). Predictors were required to 
start these computations in a prescribed way, to use a prescribed 
set of free-stream conditions, to plot the results on a standard 
form, and to report all free -parameter adjustments. Most of the 
prediction programs were set up by graduate students for operation 
on the Stanford computer, so that by the time of the conference 
considerable experience with the various methods had been devel- 
oped by the host group. This was very useful in preparing a 
paper on the morphology of the methods (Reynolds in TBLPC-1). 


* Executive Committee : M. Morkovin, G. Sovran, D. Coles. Host 
Committee : S. J. Kline, E. Hirst, W. C. Reynolds. Advisory Board : 
F. H. Clauser, H. W. Emmons, H. W. Liepmann, J. C. Rotta, I. Tani. 
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The predictors sent their computations to Stanford shortly 
before the meeting, and we compiled them for review at the con- 
ference. Comparisons were limited to three integral parameters 
of the mean flow available from all computations* the momentum 
thickness 0 , the shape factor H , and the friction factor C f . 
Mean profile comparisons were made by many authors and a few even 
made comparison with turbulence data. 

Fig. 1 shows the common comparison for the easiest flow, a 
flat plate boundary layer. On shifted scales we show H , 

C f (CF) and R 0 (RTH) vs. x for the 29 methods examined at 
the conference. The letters on the left identify the method. 

Note that all but one method is able to handle this flow reason- 
ably adequately (see TBLPC-1 for method code key). 

Fig. 2 shows the comparisons for a more difficult flow, an 
"adverse" pressure gradient (decelerating free-stream) flow. 

Note that some methods do reasonably well while others do quite 

poorly. One predictor exemplified the integrity of the confer- 
ence by producing a calculation that failed to fit his own data.' 

A small committee headed by H. Emmons studied the results 
and attempted to rank the methods. Fig. 3 shows a comparison of 
the rankings of two evaluation committee members. Methods based 
on partial differential equations are shown as P , those includ- 
ing a turbulence differential equation are indicated by P + , and 
integral methods are shown as I . The committee noted that 
several different kinds of methods performed quite well, and that 
certain methods were consistently poor. They went on to recom- 
mend abandoning the poor methods in view of the success of the 
better ones. 

While there were a number of successful and attractive 
integral methods tested at the conference, one had to be impressed 


4 



with the generality and speed of computations based on the partial 
differential equations. These schemes can be extended to new 
situations much more readily than integral methods. While inte- 
gral methods are indeed useful in certain special cases, there is 
a definite interest in use of partial differential equation 
schemes. In view of the advantages of partial differential 
equation methods, we have omitted integral methods from further 
consideration in this review. However, the development of 
adequate partial differential equations may well stimulate devel- 
opment of new integral methods based upon these equations. Two 
of the better TBLPC integral methods were developed in this 
spirit (Me Donald and Camarata, Hirst and Reynolds in TBLPC-1). 
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3. Status of Closure Experiments 


In time averaging the Navier-Stokes equations to render the 
turbulent flow problem tractable, information is lost to the 
point that the resulting equations are not closed. Additional 
equations may be derived by manipulations with the Navier-Stokes 
equations before the averaging process, but the number of 
independent unknowns increases more rapidly than the number of 
equations, and rigorous closure is just impossible. Apart from 
direct numerical solution of the unaveraged equations, about which 
we will comment later, the only hope lies in replacing some of 
the unknowns in the equations by terms involving other unknowns 
to bring the number of unknowns down to the number of describing 
equations. Such assumed relations are called "closure assump- 
tions" . 

In turbulent shear flows it has seemed most convenient to 
work with the velocities as independent field variables (rather 
than with their fourier transforms as is often done in isotropic 
turbulence). The simplest closure involves only the mean momentum 
equations, which contain unknown turbulent stresses for which a 
closure assumption must be made. This is the MVF closure (Mean 
Velocity Field). Models of this sort have been applied to a wide 
variety of flows, and work quite well for mast boundary layer 
flows of the Stanford conference. MVF methods are denoted by P 
in Fig. 3. • 

The next formal level of closure is at the level of the 
dynamical equations for the turbulent stresses, which we shall 
call mean Reynolds stress closures (MRS). There have only been 
a few experimental calculations at this level, and such closures 
are not yet tools for practical analysis. 
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An intermediate closure level using the dynamical equation 
for the mean turbulent kinetic energy (MTE closure) has dominated 
more recent calculations, and has developed to the point of 
utility as an engineering tool. MTE closures are denoted by P + 
in Fig. 3. Since MTE closures permit calculation of at least one 
feature of the turbulence, such methods work better than MVF 
closures in problems where the turbulence behavior lags behind 
sudden changes in conditions. In addition, they give more use- 
ful information for only a little additional effort, indeed for 
considerably less effort than MRS closures. They do not give 
adequate detail on the turbulent structure and do not work well 
when the structure (but not the energy) depends explicitly upon 
some effect, such as rotation. MRS closures will be needed for 
these problems (although MRS closures have not yet’ been tested 
in cases where they are really needed to obtain accurate mean 
velocity predictions). It would appear, then, that MTE closures 
will remain important for some time, serving both as useful 
engineering tools and as guides to the development of more com- 
plex models. 

Another approach that has promise for study of turbulence 
structure is the fluctuating velocity field (FVF) closure used by 
Deardorff (1970). Using the analog of a MVF closure for turbulent 
motions having scales smaller than his computational mesh, 
Deardorff carried out a three-dimensional unsteady solution of 
Navier-Stokes equations, thereby calculating the structure of the 
larger scale eddy motions. It is likely that such calculations 
will remain beyond the reach of most for some time to come, and 
that we will have to content ourselves with whatever we can 
extract from the computationally much simpler MVF, MTE, and MRS 
closures. However, results like Deardorff' s should serve as 
guides for framing closure models. 
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Truly fresh approaches to turbulence have not been frequent, 
and this review would not be complete without the mention of two 
that show promise for future research. The first is Busse’s 
(1970) and Howard's (1963.) work in fixing bounds- on the overall 
transport behavior of turbulent flows without any closure approx - 
imations . The second is the use of multipoint velocity probabil- 
ity densities (Lundgren 1967, Fox unpublished) with closure 
assumptions being made on the probability densities rather than 
on velocity moments. Neither of these schemes is presently 
developed as general analytical tool, but either could spark a 
major revolution in turbulence theory. 

In the sections that follow we will outline the theoretical 
framework of the MVP, MTE, and MRS closures, and present examples 
and commentary on applications of each method. Readers unfamiliar 
with the differential equations should consult Hinze (1959) or 
Townsend { 1956 ), In several instances I have taken the liberty 
to reformulate the constitutive models in an effort to extend 
their generality. 

Following up on the concerns expressed about invariance at 
TBLPC, I have made extensions that put the basic equations in a 
properly invariant manner. One must not read too much into this, 
however. Bradshaw (1970) cites Russell's (1961) wisdom: "A 

philosophy which is not self-consistent cannot be wholly true, 
but a philosophy which is self-consistent can very well be wholly 
false ... There is no reason to suspect that a self-consistent 
system contains more truth. " 



4. MVF Closure Theory 


The equations for the mean velocity field U. and pressure 
P in an incompressible fluid with constant density and viscosity 
are 


du. 

X 


Su ± + „ s "i..up . v a2u i . 

3t U j dx . p 3x7 3x73x7 6x. 

o X <} J J 


(4.1b) 


where R. . = u.u . . We will loosely call R. . the Reynolds 

XJ X J XJ 

stress tensor (actually -pR. . is the stress tensor). The over 

i J 

bar denotes a suitable average, and u^ is the instantaneous 
fluctuation field. 


One obtains closure through assumptions that relate the 
Reynolds stresses R. . to properties of the mean velocity field 

X J 

Ih . The most productive approach has been to use a consitutive 
equation involving a turbulence length scale, usually called the 
"mixing length". A generalization of the usual assumption is 




2 \J~2S r 

V mn mn 


lZ s u 


(4.2a) 


,30. 3u . 

where + 3 ^) 

.11 


(4.2b) 


2 

is the strain-rate tensor, q = R^ = u/u^ , and i is a 
turbulence length scale. Throughout we shall denote such 
length scales by l , often subscripted. For the special case 
of simple shearing motion, where 
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s, , = 

a. j 


(4.2) gives 


n 1 dU 

0 5 1y 


1 dU 

2 dy 


(^. 3 ) 


R. . = 
1J 


2/ _^2 dU dU 

dy dy 


,2 IdU dU 
dy dy 


q 2 /3 


q 2 /3i 


(4.4) 


Now, if the spatial distribution of l is assumed, (4.1) and 

(4.2) form a closed system of equations for the variables U. 

2 2 
and P + pq /3 . Note that the combination of pq/3 with P means 

2 

that q need not be evaluated. 


Another closure approach used at this level is generalized 


E ij - f- B i 3 - 2 Vij 


(4.5) 


where v^, is the turbulent or eddy (kinematic) viscosity. An 
assumption of the ,spatial distribution of v^, also suffices for 
closure. Occasionally these approaches are mixed. Comparison of 
(4.2) and (4.5) gives 


v_ = r \/2s s 

T V mn mn 


(4.6) 
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and consequently assumptions about & are often used to determine 
Vrj, , or vice versa. 

Mellor and Herring (1970b, hereafter referred to as MH), 
observing that (4.2) or (4.5) imply that the Reynolds stresses 
deviations from q /3 6^ are proportional to the strain rates 
(and hence that the principal axes of the stress deviation and 
strain-rate are aligned ), call these closures "Newtonian". Accord- 
ingly, we denote them by MVFN. The success of the Newtonian 
model is remarkable, especially since for even the weakest of 
turbulent shear flows the principal axes are not aligned 
(Champagne, Harris, and Corrsin 1970). 

In MVFN calculations the mixing length i is assumed in 
terms of the geometry of the flow. In a thin free shear layer, 
such as a Jet or wake, the assumption that b is proportional 
to the local width of the layer seems to work quite well, with 
something like 

l = 0.16 (4.7) 

This behavior is also used in the outer region of a turbulent 
boundary layer. Near a wall l is experimentally found to be 
proportional to the distance from the wall, and 

i = xy k is? 0.41 (4.8) 

seems to hold for smooth walls, rough walls, with modest com- 
pressibility, with transpiration, and in Just about any axial 
pressure field. 

In the viscous region immediately adjacent to a wall the 
calculations are improved if l is reduced, with 
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I = Ky[l-exp(- ^|^ -)] (4.9) 

where u* =fe w /p is the friction velocity based on the local 

wall shearing stress ? , and A + is a parameter character- 

w 

izing the thickness of the viscous region on the familiar y+=yu*/v 
scale. A + is known to depend upon both the streamwise pressure 
gradient and the transpiration velocity (for suction or blowing). 
Physical models of the wall layer can be used to suggest (4.9). 

Kays and his associates (private communication) have cor- 
related their turbulent boundary layer data to produce the A + 
correlation shown in Fig. 4 . There P Q + is the streamwise 
pressure gradient parameter 

p + = dP 

° pu*3 dx 

and V o + is the transpiration parameter 

C - V u * 

where V is the injection velocity normal to the porous wall, 
o 

Kays also modifies (4.9) by using the local shear stress ^(y) 
rather than "T w ' in u* . 


(4.10a) 


(4.10b) 


In boundary layer calculations, most workers simply use 
zonal models, with (4.9) in the inner region (which becomes (4.8) 
further from the wall) and something like (4.7) in the outer 
portion of the flow. Byrne and Hatton (1970) use a three-layer 
model as the basis for v,p assumptions. MH have used concepts 
from the theory of matched asymptotic expansions to obtain 
composite representation for i valid across an entire turbulent 
boundary layer. A typical distribution of l in a boundary 
layer is shown in Fig. 5. 
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For steady, two-dimensional incompressible boundary layers 
the MVFN equations reduce to (th = (U,V,W), x i = (x,y, z) ) 


3U , dV _ n 
33E + 37 " 0 




c>p 


oo . o / — \ 

+ (-uv) 


p cTx" 37 


and 


-uv = i L 


dU 

dy 


dU 

37 


or 


-uv 


dU 


T 37 


(4.11a) 

(4.11b) 

(4.11c) 
(4. lid) 


These equations are of parabolic type, and may be solved by a 
forward marching technique. The upstream profile U(x ,y) 
must be specified, and. the free-stream pressure distribution 
P oo (x) must be known. V(x Q ,y) is then determined by (4.11a). 

The numerical problems are straightforward but not a trivial 
aspect of a successful method. Implicit schemes have been most 
successful, although explicit marching methods can be used if the 
wall region is treated separately. 


In order to handle the rapid variations near a wall, one 
must either use a fine computational mesh in this region or else 
employ a special treatment. The variation in shear stress is, to 
a first approximation, small across this region, and the 
"law of the wall" is known to be followed by the mean velocity 
profile very near the wall for most turbulent boundary layers. 

One simple approach is therefore to patch the numerical solution 
at the first computation point away from the wall to the empirical 
wall law. 



(4.12) 
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This sets the value of U in terms of the wall shear stress 
(taken as the shear stress at the first mesh point) and y value 
at that point, and V may be taken as zero (or V 0 ) there. 

These conditions then provide boundary conditions for. the numer- 
ical solution in the outer part of the flow, and a nearly uniform 
computational mesh in the outer region is usually feasible. 

For transpired boundary layers or strong favorable pressure 
gradients the shear stress variation in the wall region is sig- 
nificant and a better analysis is required. One approach is to 
use a solution to the governing equations obtained by assuming 
parallel flow (neglecting axial derivatives, except for pressure). 
This "Couette flow" solution is obtained by analytical or numeri- 
cal solution of ordinary differential equations, and these solu- 
tions may often be precomputed in parametric form. A semi- 
theoretical wall layer treatment of this sort is very effective in 
permitting large computational steps in the streamwise direction. 

The Couette flow analysis uses the constitutive equation as 
its basis. The total shear stress in the boundary layer is 
written as 


f=( v + Vsf 

Eq. (4.13) may be integrated and expressed in dimensionless form. 


U 

u* 



T + 

(l+v T /v) 



(4.14) 


where y + = yu#/ v and T + -T/V- Thus, to develop the inner 
region solutions one needs to know the shear stress distribution 
T (y) . In the Couette flow approximation the convective terms 
are deleted, and the shear stress emerges from the momentum 
equation as 

T + - 1 + U + V 0 + + P + y + (4.15) 
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Loyd, Moffat, and Kays (1970) have found this inadequate for 
strongly accelerated flows beyond y + = 5 . Since the patching 

will take place at a much larger value of y + (perhaps around 30 - 
50), a better shear stress distribution is needed. Loyd et al 
noted that for fully asymptotic flow, where U/U^ = f(y/ 6 ) 
throughout the entire layer (such flows can be realized with 
strong acceleration), the shear stress distribution is 

y 2 

t+ =. 1 + u + v o + + P + y + [1 - J [ (^) dy] (4.16) 

and they use this expression to obtain a better shear stress 
distribution for use in the Couette analysis. These integra- 
tions are carried out at each streamwise step in the computa- 
tion to patch the inner and outer solutions. 


Recently Kays has found that improvements in the prediction 
of flows with sudden changes in wall conditions are possible if 
empirical "lag equations" are used for the parameters F* and 

+ + 1 

V Q used to determine A from the correlation of Fig. 4. 

Loyd, Moffat, and Kays (1970) use 


+ 4. + 

dP Q + P -P. 
e _ e 

dx + C x 


dV 


oe 


dx + 



(4. 17 a, b) 


with and^ Cg of approximately 3000. Here P + and V Q + 

are the actual values, and P e + and V 0 + are the "effective" 
values used in reading A + from Fig. 4, and x = xu*/ v .- 


Fine wall mesh schemes have been used to avoid this patching 
process. It is critical to use a good implicit difference scheme 
in this case. Mellor (I 967 ) developed a good linearized iteration 
technique which has since been adopted by others. 
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The approach to calculation of the temperature field and 
heat transfer follows closely the hydrodynamic calculation out- 
lined above. For incompressible flow of a fluid with constant 
and uniform properties, neglecting the input to the thermal 
field by viscous dissipation, the thermal energy equation (ob- 
tained by a combination of the energy and momentum equations) is 


de 


+ u 


be 


a d 2 9 



(4.18) 


Here 0 denotes the mean temperature and 9 the local tempera- 
ture fluctuation. The terms u .© represent transports of 

J 

internal energy by turbulent motions, and it is these terms that 
bring the closure problem. 


The common approach to the thermal problem is to assume 



Ctn 


80 


(4.19) 


where is the "turbulent diffusivity for heat", analogus to 

v,p . With knowledge of , and with from solution of the 

hydrodynamic problem, the thermal problem is closed. It is 
usually assumed that 


v I /a T = Pr T (4.20) 

• / 

where Pr,p is a turbulent Prandtl number. For gases Pr^ is 
experimentally found to be approximately 0.7-1 in typical boun- 
dary layer flows, and a constant value often suffices. More 
elaborate correlations of Pr^, with other properties of the 
flow have also been proposed (Simpson, Whitten, and Moffat 1970, 
Cebeci 1970b). The choice of Pr,p is particularly important for 
liquid metal heat transfer. 
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In examining the nature of a T and v T in the viscous 
region of boundary layers, use has often been made of an unsteady 
two-dimensional parallel flow Stokes model (Cebeci 1970b). While 
such analysis may well yield the relevant dimensionless groupings, 
and possibly a fairly reasonable form for the a T and v T dis- 
tributions, failure to consider the now well established strong 
three-dimensional unsteady features of the laminar sublayer 
(Kline, et al 1967) would seem to render quantitative results 
questionable. Since the heat transfer rate in boundary layers is 
strongly dependent on the assumptions made in this region, it 
would seem that at present the best results will be obtained 
with models having high empirical content, such as the A + cor- 
relation of Fig. 4 and the Pr T correlations of Simpson et al 
(1970). New theories based on more accurate models of the wall 
layer will probably get considerable attention. 


Though the concept of a turbulent viscosity has been dis- 
pleasing to many, one cannot deny the success that its users have 
enjoyed. An interesting interpretation of v T is obtained by 
multiplying (4.5) by S ij . , viz. 



v 


-R, .S. . 
ij 

2vS S 
mn mn 


(4.21) 


The numerator is the rate of production of turbulence energy, 
and the denomenator is the rate of dissipation of mechanical 
energy by the mean field. 
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5. MVF Calculation Examples 

Many examples of MVF calculations have now been published, 
and we shall now look at a small but representative collection. 
Readers should see the original papers for description of the 
details. 

Most published computations have dealt with boundary layers. 
The numerical techniques employed have varied considerably, and 
hence computational costs initially varied widely between pro- 
grams. But now most workers have adopted implicit difference 
schemes, with special wall region treatment as outlined above, 
and/or a linearized iteration technique (Mellor 1967 ), so that 
run times are now reasonably uniform. A typical two-dimensional 
compressible boundary layer can now be treated in voider one 
minute on a typical large computer. 

Among the pioneers and current advocates of the MVFN equations 
were A. M. 0. Smith, and his colleagues, chiefly T. Cebeci. They 
elect to specify the eddy viscosity distribution, using a form 
derived from the mixing length model in the inner region and a 
uniform value reduced by multiplication by an intermittency 
factor in the outer region. The curves marked CS on Figs. 1 
and 2 are by their method. Cebeci et. aL (1969* 1970a) have 
extended their method -to include heat transfer and compressibility. 

D. B. Spalding has been an active explorer of turbulent 
boundary layer computational methods. His early work with 
Patankar ( 1967 ) was based on the MVFN equations with mixing length 
specifications, and their complete program descriptions served as 
the seed for numerous computational efforts elsewhere. Fig. 6 
shows their computation of a wall jet flow as presented in TBLPC-1. 
This computation was among the few "more difficult" flows volun- 
tarily presented by predictors to illustrate the range of their 
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method. Spalding has now essentially abandoned this method in 
favor of MTE models. 

G. Mellor and his coworkers have used MVF closures for a 
variety of problems, and their unpublished work on the theoreti- 
cal foundations of the theory has been both educational and use- 
ful in writing this review. Mellor and Herring startled TBLPC 
by presenting two methods, one based on MVF closure and a second 
based on MTE closure; except in one case the H , 6 , and C f 
predictions by the two methods were absolutely indistinguishable, 
both being judged among the best at the conference (shown as MH 
on Figs. 1 and 2 ). Mellor (1967) has also used a MVF method to 
study certain classes of three-dimensional boundary layers, and 
Herring and Mellor (1968) have extended the method' to compressible 
boundary layers. 

Since TBLPC, interest in the MVFN prediction methods has 
spread. F. Dvorak (private communication) has been looking at 
applications to more difficult flows of interest in aircraft 
design, and has kindly provided Fig. 7 as an illustration of his 
work. With some adjustment of the eddy viscosity prescription, 
Dvorak is able to predict the growth of a boundary layer with 
tangential injection upstream and a strong adverse pressure 
gradient. This flow has two overlaid mixing layers, which 
suggests the variation in v^, used by Dvorak, though it would 
seem difficult to make really accurate calculations if the down- 
stream data were not available to guide the v T tailoring. 

The MVFN equations have been used in the calculation of 
three-dimensional boundary layer flows by Mellor (1967) and 
currently by Wheeler and Johnston (unpublished). We remark that 
the MVFN model assumes that the shear stress is aligned with the 
strain rate. In spite of the strong experimental evidence 
(Johnston 1970 ) that this does not hold, the MVFN equations work 
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remarkably well in predicting the mean velocity field in three- 
dimensional boundary layer flows where the pressure field (rather 
than the turbulent stress field) has the primary influence on the 
three-dimensionality (most boundary layers of engineering interest 
may be of this type). Fig. 8 includes integral parameter pre- 
dictions using Mellor and Herring's MVFN method by Wheeler and 
Johnston (private communication) of the flow along the symmetry 
plane in a boundary layer approaching an obstacle. Except very 
near the separation point, results are excellent. The MTE 

predictions on Fig. 8 will be discussed in Section 7. 

/ 

The prediction of turbulent boundary layer separation by 
MVF methods has not been very successful. Indeed, it may be 
appropriate to identify turbulent separation in terms of the 
turbulence near the wall, and this will require use of a more 
sophisticated model (MTE or MRS), quite possibly in their full 
(rather than boundary layer) form. 

MVFN methods have been used with some success in compressible 
flows. Fig. 9 shows a prediction of Herring and Mellor ( 1968 ) of 
the Mach number correction to the skin friction factor for a flat 
plate boundary layer. Fig. 10 shows their prediction for the 
boundary layer on a waisted body of revolution. Note that, while 
the momentum thickness is quite accurately predicted, the 
velocity profile details are in considerable error. Indeed, MVFN 
methods are often much better in predicting integral properties 
of the flow than in predicting local details. Geometrical effects 
neglected in the analysis are the probable cause of much of the 
discrepancy. 

Fig. 11 shows a prediction by Healzer and Kays (private 
communication) of the heat transfer coefficient (based on enthalpy 
difference) in an adiabatic rocket nozzle boundary layer flow, 
made with an extended MVFN method (no chemical reactions 
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considered). The accuracy of this prediction attests to the 
value of such methods in contemporary engineering analysis. 

MVFN methods have been used in contained . and recirculating 
flows, where the boundary layer approximations no longer apply. 
Spalding and his coworkers have led these efforts (see Gosman, 
et al 1969 ). The numerical treatment is critical here, for the 
equation system is elliptic rather than parabolic, and the entire field 
must be solved simultaneously. Computational times are consequently 
considerably longer, with several minutes being required for a 
typical flow. Recently Chin and Seban (private communication) 
studied an improvement of Spalding's upwind difference treatment 
as applied to the flow in a cavity under a turbulent shear flow. 

The results of their computation are shown in Fig.' 12. They 
used a simple wall region patching treatment, with a linear mixing 
length near the walls, a uniform mixing length in the central 
region of the cavity, and a constant mixing length in the external 
shear layer. The computational mesh was 41 x 41 in the cavity, 
with closer spacing near the walls. In order to obtain convergence 
in the solution of the difference equations, over 1000 relaxation 
iterations were required, and the computation took 20 minutes on 
a CDC 6400 computer. While the velocity distribution in the 
central cavity is predicted very well, the heat transfer from the 
cavity bottom is not. Seban (private communication) states that 
an improved wall region treatment is required, but that the relax- 
ation iteration became nonconvergent when this was tried. He 
suggested that perhaps the time -dependent MVFN equations would 
have to be solved in order to compute the final steady-state flow. 

MVFN equations have not been tested in Aery many time- 
dependent flows, for there is practically no comparison data. 
Moreover, the computation costs skyrocket with every added dimen- 
sion. However, if the time -dependence is periodic, a fourier 
analysis can be used to reduce the problem to a sequence of steady 
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problems. If the flow is parallel and the periodic component 
takes the form of streamwise travelling waves of small amplitude, 
then the MVFN equations may be reduced to ordinary differential 
equations for the periodic disturbance similar in structure to 
those used in analysis of the stability of laminar flows. We 
have been looking at the results of such computations for periodic 
disturbances in shear flows and for flows over waving boundaries. 
Our experimental observations of small periodic disturbances in 
turbulent channel flow (Hussain and Reynolds 1970) indicate that 
a dispersion relation exists between the frequency and stream- 
wise wavenumber of disturbance eigenfunctions. Pig. 13 shows our 
predictions for this relationship as compared with our experi- 
mental data. The predictions were made using the eddy viscosity 
distribution calculated from the mean velocity profile, a fine 
wall mesh, and (4.5) in the time -dependent MVFN equations. Note 
that the MVFN model seems to work well in this unsteady flow. 

We have also applied this approach to flows over waving 
boundaries, and in particular to Kendall's (1970) flow and 
Stewart's (1970) flow. In neither case did our predictions 
agree with the measurements; Davis (1970) used a similar MVF 
model with curvilinear coordinates, apparently with greater 
success. An experiment on turbulent channel flow with a waving 
wall has just been completed in our laboratory. The wave-induced 
wall pressure oscillation is predicted fairly well by MVF theory for 
upstream-running waves, but not at all well for downstream-running 
waves. This suggests that the MVF model is weakest in flows with 
a "critical layer", i.e. a point where the mean velocity matches 
the wave speed. The ability to predict such flows by MVF 
methods would seem questionable, in view of the rapid changes in 
strain rate to which the turbulence is subjected, and in all 
probability a MTE or MRS method would work much better. We 
intend to explore calculations along these lines. 
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MVF methods fail in any flow where the nature of the 
turbulence is altered by some parametric effect, such as rotation, 
which does not appear parametrically in the equations of mean 
motion. Such effects can be included in MVF methods only by 
alteration of the l or v,p specification, and hence MTE or 
MRS methods are clearly to be preferred for such cases. 

The most ambitious application of MVFN equations has been 
to atmospheric general circulations. The National Center for 
Atmospheric Research (NCAR) has developed an elaborate model in 
which the velocity components, temperature, and humidity are 
calculated over the entire earth. The goal is to obtain an 
accurate 14-30 day weather forecast. The computational mesh 
involves six vertical layers and 5° grid spacing at the equator, 
with fewer points near the poles. The horizontal grid therefore 
varies from about 500 km to 100 km on a side. A turbulent 
viscosity model is used to handle sub-grid scale turbulence. The 
effects of sun, show, water, mountains, and precipitation are 
simulated. The main features of global weather patterns are 
reproduced. The dearth of field data make quantitative comparison 
difficult, and initialization almost impossible. Kasahara (I 969 ) 
reports that "better" results are obtained with a 2.5° mesh, but 
with present computers a 24-hour computation requires about 24 
hours with this finer mesh. Fig. 14 shows an NCAR computation 
from the 5° model. It seems quite possible that such calculations 

will someday become a routine part of our weekly weather forecasts, 

\ 

though refinements in the physical model may be required. 
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6. MTE Closure Theory 


The MVF equations assume that the turbulence adjusts immed- 
iately to changes in mean conditions, and that a. universal re- 
lationship exists between the turbulent stresses and the mean 
strain rates. To avoid these assumptions one must include dif- 
ferential equations for the Reynolds stresses (called "dynamical" 
or "transport" equations). MRS closures use these equations; 

MTE closures are somewhat simpler, and employ a single equation 
for the turbulent kinetic energy in conjunction with some 
constitutive or structural equations relating the turbulent 
stresses to the turbulent kinetic energy. Thus, MTE methods can 
to a degree handle the delayed response of turbulence structure 
to sudden changes in mean conditions, and are now being studied 
by several groups for use in such problems. 


Equations for the R. . may be developed from the Navier- 

X J 

Stokes equations (Townsend 1956; Hinze 1959). These are 


dR. . dR. . dU dU 

~St + u k = “ R ik ' R jk 3x£ 


" ^ (u i u j u k ) " p [ s!t ( u i p) + (u j p)] 


_ du. du. 

+ £ (-5-X + —Li) + V. . 

p v dx . dx . 1 ij 
J x 


( 6 . 1 ) 


Here V. . is the viscous term to be discussed shortly. A con- 
J 

traction of these equations gives the equation for the turbulent 
kinetic energy. With q = u.u. , this may be written as 
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( 6 . 2 ) 


^6^ + \ ' ' R nc 33T ■ sir K u iV 2 ) 


,x k 


P st^ (u k p > + v ii /2 


The first term on the right is the "turbulence production". 


The more common form of V. . is 

X J 

d 2 R,, 

^dx k ‘ 2 °^ij 


where 


Then, 


V. . = v 
i J 


(6.3a) 


A du. du. 

A, -V -s-i 3-i 

U ox k ox^ 


(6.3b) 


where 


V_. 


e 2_2 
V dx/dXv 


ii _ a g q g /2 


k k 


- id 


(6.3c ) 


A.. 

*11 


(6.3d) 


This form is appealing because the first term in V.J./2 can be 
interpreted as a "gradient diffusion" of turbulent kinetic 
energy, and the second is negative-definite (suggestive of 
"dissipation" of turbulence energy). However, the rate of entropy 
production is proportional to 


du, du . du . 
v + 3 ?) 3337 ^° 


(6.4) 


Properly € is called the "dissipation" but not id • We might 
call $0 the "isotropic dissipation". 
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A second form of is 


v ij = v 33T K °jk + Vikl 


du . du . 

- v t°jk ssr + ff ik ss? 1 


where 


du . du . 

r ij ~ dx . + dx. 


For which (6.5a), 


v u 


~2~ ~ v 5x~ [u i°ik ] ' e 


(6.5a) 


(6.5b) 


(6.5c) 


The appearance of e makes this form appealing, even though the 
first term can no longer be interpreted as "gradient diffusion". 

MTE methods require closure assumptions for the last three 
terms in (6.2) and there has been heated debate on this point. 
There seems to be universal agreement that the dissipation term 
should be modeled by the constitutive equation 



where l is a "dissipation length scale", and C is a function 
of the dissipation Reynolds number R € = q ,/ / v , C being con- 
stant for R_ » 1 . Note that for R_ » 1 e is independent 
of v ; this is a reflection of the belief that the small scale 
eddies responsible for the final dissipation of mechanical energy 
can handle all the energy that is fed to them from and by larger 
scale motions, and hence the larger eddies control the dissipation 
rate. The spectral transfer process for R £ > 1 results from the 
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inertial non-linearity* which suggests (6.6). The remainder of 
the viscous term is only important very near a wall; though it 
is strictly incorrect* reasonable results have been obtained by 
taking < 0 = * * and writing 


V ii 


/ 2 = 


'’if#-* 


MH's model is more complicated (see 8.10* 6.22), 


(6.7) 


The main MTE argument stems over the treatment of the 
pressure-velocity correlation term and the triple velocity 
correlation term. One widely used approach is the "gradient 
diffusion" model* where one sets 


( u k p + u k 


u i u 


-) :« 


Q V T 


Bq 2 _/2 


( 6 . 8 ) 


where is a constant (or specified function). There is 

strong feeling in some quarters that this model ignores the 
dominance of transport processes by large scale eddy motions. A 
generalization of a "large eddy transport" model (Bradshaw* 
Ferriss* and Atwell 1967) is 

. « m 

< u kp + u k ^r) = Gq\ (6.9) 


where G is a constant (or specified function) and is a 

global vector velocity scale characteristic of the large eddy 
motions. The choice for this closure is of considerable importance; 
neglecting the viscous diffusion terms* the equation system based 
on (6.8) is of elliptic type* while with ( 6 . 9 ) the system is 
hyperbolic. This mathematical difference is suggestive of sub- 
stantial physical differences in the model. Both approaches 
have been used quite successfully, however* and it. is not easy to 
make a strong case for either solely by testing against experiments. 


27 



2 

Having closed the q equations, one must relate the Reynolds 

P 

stresses to q in order to have a closed system. Again two 
approaches have been hotly debated. The more common approach 
uses the constitutive equation (4.5), together with an additional 
constitutive equation relating the turbulent viscosity to the 
turbulent kinetic energy, 

v T = qiF(R T ) (6.10) 

Here l is a turbulence length scale, and F describes the 
dependence upon the turbulence Reynolds number R^ = q i/v 
with F = const, for R^ » 1. The length scales l and 
must be specified (either algebraically or through a differential 
equation) to close the equation system. Since use of (6.10) 
again implies Newtonian behavior, we shall refer to this MTE 
closure as MTEN. 

Observing that the Newtonian structure is never observed in 

turbulent shear flows, but that persistently strained flows 

apparently develop an "equilibrium structure", Bradshaw prefers 

2 

to relate the Reynolds stresses directly to q . A generaliza- 
tion of his constitutive equation is 

R ij - ai/ ( 6 -U) 

case of 
Townsend 


( 6 . 12 ) 


where a. . ' depends upon the type of strain. For the 
**• 

pure shear (4.3), a reasonable form of (6.11) is (see 
1956, Champagne, Harris and Corrsin 1970) 

( 0.48 -0.16 0 ' 

-0.16 0.26 0 

0 0 0.26 , 


28 



Lighthill (1952) suggested a general form which gives a l2 = -0.16 
but does not correctly represent the diagonal terms. 


a. . = | 6. . 

ij 3 iJ 


°_± S i,1 


mn S mn 


(6.13) 


We will denote MTE closures involving an assumed turbulence 
structure (e.g. 6.12) by MTES 


One would like to assume constant values for a^ in thin 

shear layers. However, on a symmetry axis i a pipe or free jet 

flow, where R 12 = 0 , a l2 = 0 , and hence to use (6.11) in such 

flows one must specify a variation in a^. . Hence, one must 

have a good "feel" for the flow to obtain a good prediction. 

This requirement for intuition is less important in simpler 

boundary layer flows, where a uniform value of a. . produces 

J 

reasonable results. 


There is a more fundamental objection to the MTES idea. 
Recently Lumley (1970) has argued that the homogeneous flows 
upon which (6.11) and (6.12) are based do not really reach 
equilibrium, and that instead the turbulence time (and length) 
scales continually increase. Champagne, Harris, and Corrsinfe( 1970) 
experiments confirm this expectation. Hence, a structural model 
cannot be fully correct in homogeneous flows. 

MTEN and MTES closures both fail in the case of a sudden 
removal of the mean strain rate, where it is known that a very 
slow relaxation of the structure towards isotropy takes place. 

The MTEN model instantly becomes isotropic, which the MTES model 
retains a permanent structure (unless one twiddles with the a. . ). 
This may not be a serious objection as long as these methods are 
used in shear flows having reasonably persistent strain. 
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To summarize, the MTE closures commonly employed use one of 
the following two forms: Using (4.5), (6.6), (6.7), (6.8), and 
(6.10), with £ = i € , 


+ u j Sv 


T S i3 S id 




+ 




(6.14a) 


v T = qiF 

Or, using (6.6), (6.7), (6.9), and (6.11) 



(6. 14b) 


(6.15) 


For (6. 14) values or distributions for C , F , and Nq must be 
assumed, while for (6.15) values or distributions for a.. , C , 

G , and • are needed. Both forms require an assumption for the 
spatial distribution of the length scale £ , The terms with v are 
not important except very near walls, and are often neglected in 
the outer flow. 


Most computations have used length scale distributions of 
the sort described in Section 4 above. Recently there has been 
some interest in using a differential equation for l , and the 
most extensive test of this approach has been by Rodi and Spalding 
(1970), and Ng and Spalding (1970). Their length scale equation, 
which is based on a spectral transport equation (Rotta 1951), can 
be generalized with slight modification as 
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°iV s ij s i 3 


c 2^ 3 + c 4 sIt^t sl:^ 2 / 2 )] 

J J 


(6.16) 


Spalding and his coworkers are able to obtain very good predic- 
tions of a variety of boundary layer and free shear flows using 
(essentially) (6.16) to determine l , provided some adjustments 
in C 2 are made near solid walls. 

Gawain and Pritchard (1970) proposed a more complicated 
hueristic integro-differential equation for turbulent length 
scales. In effect their local length scale is determined by the 
mean velocity field in the region of the local point. The two- 
point tensor 


£ ) = ^-(x + 6)uj(x - %) (6.17) 


is used to define the length scale, 

i 2 = JR ±± £ 2 dV//ft i:L dV (6.18) 

where d-V demotes a volume integration. A form for is in 

effect assumed in terms of the mean velocity field,- and the inte- 
grations are performed to obtain l . This length scale is then 
used in a MTEN calculation method, where reasonably accurate 
results are reported for plane Poise.uille flow and for an axi- 
symmetric jet flow. 

Harlow and Nakayama (1969)* noting that the length scale 
will be used to determine the dissipation, proposed a closure 
model for the exact differential equation for $8 derivable from 
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the Navier-Stokes equations. They experimented with the use of 
this equation in MTEN closures. The Los Alamos group (private 
communication) has now abandoned the MTE closure in favor of MRS 
closures, which also use the c0 equation for inference of length 
scales (see 8.13). They refer to the {& equation as a "dissipa- 
tion" equation, which as we have noted is not strictly correct. 

Hanjalic, Jones and Launder (1970) have used a dissipation 
model equation to study a variety of boundary layer flows in an 
extended MTEN model. Their formulation is purported to work in 
the viscous region, eliminating the need for wall-solution patch- 
ing ( see 8.14). 

The interest in and activity with dynamical equations for 
the dissipation (or length scale) suggests that such equations 
will shortly become an important and well advertised feature of 
MTEN prediction methods, and probably of MRS methods as well. 

The dissipation equation is discussed in greater detail in Section 

8 . 


The boundary layer form of (6.I5) is (neglecting viscous 
terms for the outer region) 

(|t + u j S77) q2/2 = aq2 " cR I " ^ Gq2Q 2^ ( 6,19a ) 

J 

a ■ = Ia l2 (6.19b) 

'Z'/p = - uv = aq 2 (6.19c) 

Bradshaw, Ferris, and Atwell (1967) use ( 6 . 19 ) to derive a dif- 
ferential equation for the turbulent shear^ st ress ^ . The 
transport velocity Qg is taken as \J ^" max / P * where ^ ax is 
the maximum value of T(y) in the boundary layer. G and &' are 
prescribed as functions of the position across the boundary layer. 
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and a is essentially taken as constant. Together with (4.l0a,b) 
(6.15) gives a closed set of equations for U , V , and 'C this 
system is of hyperbolic type, with three real characteristic lines 
Bradshaw, Ferris s and Atwell construct a numerical solution using 
the method of characteristic si it can also be done using small 
streamwise steps with an explicit difference scheme (Nash; 

Wheeler and Johnston private communications ) . There is a great 
physical appeal to the characteristics, especially since it is 
found that the solutions along the outward going characteristic 
dominates the total solution. This may well be connected with 
physical observations on the nature of turbulent boundary layers 
(Kline, et al 1967 ). 


The boundary layer form of (6.14) is (neglecting viscous 
terms for the outer flow) 


( 3 I + U i = V tO h + 3^Vt^T ] (6.20a) 


V T = qiF 


(6.20b) 


Then , together with (4.5) and (4.1 a,b) this gives a closed system 
of equations for U , V , and q , provided C , , and l 

are specified. This system is of parabolic type. 


Equations (6.19) and (6.20) do not hold in the viscous 
region near the wall. One must either modify these equations to 
include viscous effects, or else use special solutions as dis- 
cussed in Section 4 in this region. Experiments reveal a nearly 
uniform distribution of q in the wall region, except very close 
to the wall (yu*/v<20) . Moreover, the value of q/v* seems to 
be nearly universal, with 

q ~ 2 . 5u* a: ^ u* (6.21) 

Jx 
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This has been used as a "wall" boundary condition for the solu- 
tion of (6.15) or (6.20). 


MH prefer to use equations containing the viscous terms to 
calculate the inner region directly. Now the manner in which 
Vii is written and modeled becomes important, and MH's version 
of the MTEN equation can be written as 


(si + 


U. 


sirHf-) = 2 


V T S ij S ij 



( 6 . 22 ) 


The 5/3 factor yielded by MH's treatment of V.^ is a main point 
of the difference with others, and MH's rationale seems most 
cogent. MH then use a fine mesh near the wall, with l and l 
varying linearly in the wall region, and being uniform in the 
outer flow. Some of their predictions are discussed in Section 7. 

MH also examine MRS closures, and show how the MTEN closure 
results- from the MRS equations with the additional assumption of 
small departures from isotropy. While this approach is academi- 
cally interesting, even the most weakly strained flows are far 
from isotropic (see Champagne, Harris, and Corrsin 1970), and 
hence the main selling point for MTE methods is that they work 
very well for predicting a wide class of turbulent shear flows. 
Examples are given in the following section. 


MTE boundary layer equations require the same upstream infor- 
mation as for MVF computations, plus the upstream distribution 

and free-stream distribution of q . Normally the free-stream 
turbulence is set zero, but the effect of non-zero free-stream 
turbulence can be incorporated in a MTE calculation. If a dynami- 
cal equation for the length scale is used, then the upstream, free- 
stream, and wall boundary conditions for l must be given. The 
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upstream i distribution can be drawn using the ideas in Section 
4, At the boundary layer edge di/dy = 0 seems appropriate. In 
the wall region 1 = 0 if the calculation is carried to the wall, 
and i = Ky if the mesh computation is patched to a wall region 
solution at the innermost mesh point. The need for this turbu- 
lence information makes MTE methods somewhat more difficult to 
use, but the ability of a good MTE method to predict more severe 
test flows may make the extra effort worthwhile. 

In the so-called "log" region of turbulent boundary layers, 
the turbulence energy is essentially determined by a delicate 
balance between the production and dissipation terms in (6.2). 

With q = u*A and i = xy (see 4.8, 6.21) a balanpe between 
the first two terms on the right in ( 6 . 19 a) gives 


dU = /C x uf; 
dy '•aK-' xy 

while for (6.20) one has 



(6.23) 


(6.24) 


Hence, both models will give the proper logarithmic velocity 
profile (4.12), provided the parenthetical coefficients are 
unity in each case. 


Heat transfer predictions made using MTEN closures have 
employed the models described in Section 4. The hydrodynamic 
calculation yields , and (4.20) is then used in (4.18) to 
construct the temperature field. 
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7 . MTE Calculation Examples 


The MTES approach has been advocated by Bradshaw and his 
coworkers (Bradshaw, Ferriss, and Atwell 1967, Bradshaw et al 
19 66 et seq. ). Their predictions for TBLPC must be judged among 
the very best. The agility of Bradshaw's MTES method to predict 
severe flows was demonstrated at TBLPC by their results for the 
Tillmann ledge flow, a boundary layer flow immediately downstream 
of a turbulent reattachment point (judged the most difficult 
TBLPC flow). Fig. 15 shows their predictions at a point down- 
stream in this flow, including results for two drastically dif- 
ferent initial shear stress distributions. Note that the pre- 
dictions are very insensitive to the initial (upstream) shear 
stress distribution. In spite of the modest disparity between 
the measured and predicted velocity profiles, the predicted 
momentum thickness and skin friction were in considerably better 
agreement with the data than were the MTEN and MVF predictions. 
Indeed, one left TBLPC with the feeling that Bradshaw and Ferriss' 
MTES method was likely to be the wave of the future. 


Nash (1970) has used a combination of Bradshaw's MTES ideas 
and a Newtonian assumption to treat three-dimensional turbulent 
boundary layers. Nash takes Bradshaw's structural assumption for 
the total shear stress vector, 

^uv) 2 + (wv) 2 ' = aq 2 (7.1) 

but then uses the Newtonian approximation 


27 du/dy 
uw ciw/oz 


(7.2) 


which assumes alignment of the shear stress and strain rate 
vectors. The evidence is clear that (7.2) does not hold, yet it 
fortunately is worst in flows with strong spanwise pressure 
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gradients, where the pressure gradients and not the shear stress 
control the mean velocity field. In Pig. 8 we show Wheeler and 
Johnston's (unpublished) calculation for the boundary layer along 
a plane of symmetry approaching an obstacle. The integral para- 
meter predictions by MH's MVFN method, Nash's MTES/N method, and 
Bradshaw's MTES method are almost identical. Nash's (1970) own 
calculation for a point off of the symmetry plane in a similar 
flow is shown in Fig. 16. Mellor (1968) made a similar calcula- 
tion with his MVFN method with comparable results. 

Bradshaw (TBLPC-1) himself extended his MTES method to 
three -dimension boundary layers, using the basic ideas to propose 
model equations for the vector sum and ratios of the two primary 
stresses, -uv and -wv . Johnston (1970) has compared the 
result of predictions by this method with his own data for an 
infinite swept flow. In particular, data show that the stress 
vector does not align with the strain-rate vector as the Newtonian 
closures assume. Hopefully Bradshaw's structural model would work 
better on this flow but Johnston's calculation shows that the 
angle of the shear stress vector is predicted quite poorly, 
although the mean velocity is predicted quite well. It is unlikely 
that MTE methods will ever predict this structural difference well, 
and one might hope that MRS methods will do considerably better. 

Me Donald and his associates (unpublished) have developed a 
MTES method following the lines of their integral method (TBLPC-1), 
and are using this method in a variety of boundary layer flows. 

They are also treating boundary layers using the full equations 
in order to study boundary layers near separation. 

There has been considerable activity with MTEN computations. 
Beckwith and Buahnell (TBLPC-1) presented partial results from 
their MTEN method at TBLPC, and have since continued with its 
development. Spalding and his associates pushed ahead with MTES 
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program development. MH have added to the theoretical framework 
through their application of the method of matched expansions to 
the selection of the length scale distribution functions, and by 
showing how the MTEN equations arise as a limiting case of MRS 
equations for nearly isotropic turbulence. 

The ability of MTEN calculations to accurately predict the 
mean velocity field and turbulence kinetic energy distribution is 
demonstrated by Figs. 17 from MH's TBLPC contribution. Their use 
of a fine computational mesh near the wall is reflected in their 
accurate prediction of the inner regions. 

The MH MTEN predictions at TBLPC were among the bestj we 
again note that these predictions were identical with those of 
their MVF method for all but one TBLPC flow. Hence, for flows 
not too rapidly shocked by changes in free-stream or wall condi- 
tions, consistent MVF and MTE treatments may be expected to yield 
nearly the same results for the mean velocity and integral 
parameters. Of course, only the MTE calculation yields the 
turbulence energy distribution directly. Fig. 18 shows a MH MTEN 
calculation for a boundary layer responding to a sudden removal 
of adverse pressure gradient. Five years ago this would have 
been regarded as a "difficult" test flow, but we see that MTEN 
methods now handle it reasonably well. 

The ability of MTEN methods to handle sudden changes in 
boundary conditions is evidenced by recent (unpublished ) calcula- 
tions by Kays and his coworkers (see Loyd, Moffat, and Kays 1970). 
They have modified an early Spalding MTEN program to the point 
where it successfully predicts the heat transfer behavior of 
incompressible turbulent boundary layers with strong pressure 
gradients and with wall suction or blowing. With sudden changes 
in pressure gradient or blowing the heat transfer coefficient 
(Stanton number) changes rapidly, and such calculations are more 
difficult for MVF methods. 
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For boundary layers the pressure gradient is conveniently 
represented by the parameter 


K = 




( 7 . 3 ) 


Fig. 19 shows a prediction by Kays (unpublished) of the heat 
transfer to a boundary layer undergoing strong acceleration 
followed by a relaxation to zero pressure gradient. Note that the 
sudden jump in Stanton number as the acceleration is removed is 
predicted quite well. Fig. 20 shows another Kays calculation for 
an accelerated boundary layer, with blowing beginning midway 
through the accelerated region and continuing through the relaxa- 
tion to zero pressure gradient. Fig. 21 shows a prediction for 
an accelerated boundary layer with blowing, with transpiration 
terminated upstream of the removal of acceleration. The remark- 
able success of these calculations suggests that MTEN methods 
are now developed to the point of utility as tools for engineer- 
ing analysis. 


The MTE methods include a calculation of the turbulence 
energy, and hence one may study the effects of variable free 
stream turbulence. Kearney, et. al. (1970) have compared such 
predictions with their data, and Fig. 22 shows a typical result 
for strongly accelerated turbulent boundary layer. 

The MTEN methods have been applied to free shear flows to a 
limited degree by Spalding and his coworkers. Fig. 23 shows 
predictions by Rodi and Spalding (1970) for the asymptotic plane 
jet. This calculation was made using their model equation for 
the turbulence length scale. Gosman et al (19^9) have documented 
the Spalding MTEN program in detail, and advocate its application 
to heat and mass transfer in recirculating flows. Readers should 
be aware that such programs are under continual development, 
which should not prevent their use in engineering analysis. 
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8. MRS Closures 


In order to compute the structure of the turbulence (i.e. 
the R. . ), one must employ the dynamical equations for the R. . 

•** J 1 J 

(6.1). This has been the subject of considerable recent interest, 
though only a few computational experiments have been carried out 
and a truly "universal" general theory has yet to be established. 
We can expect considerable activity on this front, and review the 
current status here with this in mind. 


The problem is again to set up a satisfactory closure 
structure for the unknown terms in the dynamical equations, here 
the equations for R^. . Some variation in approach is already 
evident, and we may expect some interesting debate on the choices 
over the next several years. 


Examination of (6.1) shows that the R. . equations contain 

J 

a pressure-strain-rate correlation term that vanishes in the 
contraction (6.2). The effect of this term must therefore be to 
transfer energy conservatively between the three components 
R 11 * R 22 * and R 33 * and is generally believed that this 

transfer tends to produce isotropy in the turbulent motions. 
Modelings of this term should incorporate this feature. A 
plausible model of this term (Rotta 1 951 ) , supported somewhat by 
the data of Champagne, Harris, and Corrsin (1970) is. 


P. . 
ij 


8u. 

= 

P ^3 


= 


„ 2 

f <*- 


ij 


- R. .) 

1.1 ' 


(8.1) 


An objection to this model rests on the observation that the 
fluctuating pressure field is given by a Poisson equation. 


1_ 8 2 p 

p BxTcTxT 


-r 3 [ U . U . -U . U . -U . U . -U .U . 1 

ox .dx . L l ,i l j 1 j J 1 J 


( 8 . 2 ) 
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This suggests that the . model should contain terms arising 
from interactions between the mean and fluctuating velocities, 
and should somehow reflect the dependence of the pressure fluc- 
tuations on distant velocity fluctuations. Daly and Harlow 
(1970, hereafter referred to as DH) have attempted to include 
these effects in a complex closure approximation still in an 
experimental stage (see 9* 8 a). Other MRS closure calculations 
have all used (8.1). Some new suggestions are explored in 
Section 9* 


The pressure-velocity terms have been modeled in all MRS 
computations of which I am aware by extensions of the gradient 
diffusion model ( 6 . 8 ). Donaldson and Rosenbaum ( 1968 , hereafter 
referred to as DR) use 




(8.3) 


DH use a similar expression with a more complex coefficient. MH 
suggest 



(8.4) 


Various forms of gradient diffusion models have been sug- 
gested for the triple velocity term. DR use and MH accept 


u i u j u k = 


- 9 ^ ( 


as. . Sr,. Sr,, 
r-ii + ^ 


) 


(8.5) 


The DH representation may be cast as 



If the objections to a gradient-diffusion approach are valid, 
one should presumably use an extension of (6.9). A possible 
large-eddy transport model is 


u i u A = R ij e k + + Vh 


(8.7) 


The viscous terms have also been handled in different ways. 
DR and DH use in the form (6.3). Following Glushko (I965), 

DR take 



( 8 . 8 ) 


DH put 



_ 2D 
- q 2 R iJ 


(8.9) 


and use another differential equation for D * D^ =#8/2v . MH, 
invoking arguments of local isotropy and using kinetic theory as 
a guide, propose using ( 6 . 5 ) with 


u . 0. 


kj 


+ u j a ki 


8r., &r.. 

.lk . lk 

3x^~ + 


dR, . 


( 8 . 10 ) 


and 


v 


^ u i 


0 . , 

ik 




( 8 . 11 ) 


In order to complete the closure, the various length scales 
in the models above must be prescribed or related to the other 
independent variables through a differential equation. DH use a 
dynamical equation for & , derived exactly from the Navier- 

Stokes equations and then closed by assumptions. The equation 
will be discussed shortly. DH are now considering the use of two 
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length scale equations for the dissipating and energy containing 
eddies . 

MH have given considerable thought to the MRS closure, and 
show how the MTEN equations emerge from MRS equations if it is 
assumed that the turbulence structure is nearly isotropic. 

Calculations with MRS closure models have been carried out 
by DR, DH, and Harlow and Romero (1969)* Harlow and Romero used 
the model with moderate success to study the distortion of iso- 
tropic turbulent (see Section 9)» DR considered plane turbulent 
boundary layer flow in zero pressure gradient, and specified 
what seem to be reasonable length scale distributions for this 
calculation. DR's prediction of the mean velocity profile is 
good (but not better than a good MVF or MTE calculation); their 
predicted turbulent stress distributions, shown in Fig. 24 are in 
substantial agreement with experiments. DH studied plane 
Poiseuille flow with a more complex model, obtaining less satis- 
factory results. The DH model is not accurate near the wall, and 
is currently undergoing further extension and adjustment. Hirt 
(I 969 ) gives a useful summary of the thinking behind developments 
in the Los Alamos group. 

Let's now consider the "dissipation" transport equation. 

The dynamical equation for $ is derived by differentiating 

the momentum equation for u. with respect to x., multiplying 

J - 

by 2vdu./ox. and averaging. The result is 
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( 8 . 12 ) 


Now, to obtain closure one must propose models of all the terms 
between the braces on the right hand side of (8.12). This 
requires a considerable amount of courage as well as insight; 
there is no direct experimental evidence about any of the terms, 
and one can really only conjecture as to their effect. Lumley 
(1970) has used some rational reasoning for the special case of 
homogeneous flows (Section 9 )- DH used some qualitative ideas 
about the effect of each term and proposed a model of (8.12), viz. 
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+ 



S b 

~ 33T 
q j 
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) 



q 2 V 



( 8 . 13 ) 


Here b^ and are "universal parameters, all with values 

near unity (or possibly equal to zero)", and F is .a function 
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of the turbulence Reynolds number. DH also in effect make the 
assumption that (d = e through their treatment of the R ±± 
equation. 


Hanjalic, Jones, and Launder (1970) propose a model of 
(8,12) which can be generalized as 


S0 +U 

3t + u j sir 


dx^ 
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(8.14) 


and are 


■ ^ cuiu ^2 ^ LUAiDoauooj anu x ^ j j. 2 

functions of the turbulence Reynolds number. They report that 
"encouraging" results are obtained when this equation is used 
in an MTEN computational scheme. Clearly the use of such equa- 
tions is presently quite experimental, and (8,l4) i s given here 
to illustrate the rather substantial differences in ideas as to 
how best to model (8.12). 


For the special case of homogeneous flows at high turbulence 
Reynolds numbers (8.I3) and (8.14) do have a common form. Both 
may be written as 


3t 


+ u. 


d£_ 

3x7 

J 



+ c. 



(8.15) 


where (P is the rate of production of turbulence energy. This 
form is probably quite adequate for homogeneous flows (Section 9). 


Lumley (1970) has studied the distortion of homogeneous 
turbulence by uniform strain using a limited MRS closure. In 
homogeneous flow the R. . equations become 
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( 8 . 16 ). 


Lumley closes by taking 


du. du . 


6 . . 

e -i*L 


(8.17) 


p ( u j IxT + u i §* 7 ) - f ' R itj ) ( 8 - 18 ) 


where T is a time scale of the turbulence (compare 8.1). He 
further assumes that the time scale is related to the dissipation 
rate by 

2 


T = 


C^q 

~2e 


(8. 19) 


which is equivalent to (6.6). The dissipation rate is in turn 
described by 


de _ 
dt 



( 8 . 20 ) 


as deduced by Lumley from scaling arguments • based on (8.12) 
(compare 8.15). 

Lumley has solved the equation system for homogeneous shear, 
and compared the results with homogeneous strain and homogeneous 
shear experiments. Lumley' s model predicts that the time scale T 
grows without bound, so that homogeneous flows can never attain 
an equilibrium structure . Champagne, Harris, and Corrsin's (1970) 
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experiments are consistent with Lumley's notion, but Lumley's 
model does not predict the observed structure very well. Some 
improvements on Lumley's model based on (8.15) are suggested in 
Section 9* 

It does seem clear that equilibrium is never obtained in 
homogeneous flows. In inhomogeneous flows the transport features 
apparently act to set the equilibrium structure. MTES methods 
really shouldn't work in homogeneous flows, and we might well be 
suspicious of methods when the "universal constants" are obtained 
by tests against such flows. 
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9. Some New Ideas for Homogeneous Flows 

/ 

It became apparent in preparing this review that too little 
attention has been given to systematic development of the 
closure model. The approach has been to construct a comprehen- 
sive model, with numerous universal constants, and then to select 
these constants by optimizing the average fit to a number of 
selected flows. A more systematic approach would be to develop 
the closure model in a step-by-step approach, working gradually 
through a heirarchy of experimental flows. 


In order to develop some feeling for what might be accom- 
plished, I examined the following approximations to homogeneous 
flow: 


(1) Decay of isotropic turbulence (Townsend 1956) 

(2) Return to isotropy in the absence of strain or 
shear (Tucker and Reynolds 1968, hereafter 
referred to as TR) 

(3) Development of structure under pure strain (TR) 

(4) Development of structure under pure shear 
(Champagne, Harris, and Corrsin 1970, hereafter 
referred to as CHC). 


The starting point of this analysis was the dynamical equations 
for the turbulent kinetic energy and the dissipation. For homo- 
geneous flows these equations are 



(9.1) 
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(9.2) 


ds _ c e 2 

™ “ ' 1 7 


+ C 


(fe 

2 ? 


Here (P = -R. .S. . is the rate of turbulence energy production. 

X J X J 

Equation (9*1) is exact, and (9.2) is the form in which the dis- 
sipation equations of DH, Hanjalic, Jones and Launder (1970), and 
the length scale equation of Ng and Spalding (1970) can be 
expressed (see 8.15 ). For very weakly strained flows, Lumley 
(1970) developed the term with = 4 from first princi- 

ples, and neglected the Cg terms in his weak strain model ( 8 . 20 ). 
Hanjalic et al suggest C^ = 4 and Cg = 3*2 . Ng and Spalding's 
empirical flow fitting is equivalent to = 3.9 and Cg = 3.3 • 
DH use C 1 = 4 and Cg = 2 . 

If one considers the decay of homogeneous isotropic turbu- 
lence with zero strain for large turbulence Reynolds numbers, and 
models the dissipation by 

e = q 3 /^ (9.3) 


(9.1) and (9.2) produce 
dq 2 /2 = _ q 3 

dt (C 1 ~3 )& 


(9.4) 


JpE - (C r 3)q (9.5) 

2-1 

Now, experiments indicate that q ^ t , i t which requires 
C. = 4 . This seems a clear choice. 

To investigate Cg I calculated the distribution of (P from 
the data of CHC, and carefully determined an initial value for e 
from the experimental q 2 distribution (taking the starting point 
at x = 5 ft in their experiments). The differential equations 
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(9.1) and (9.2) were then solved numerically for different values 
of C 2 ; Cg =2 is clearly preferred (Fig. 25a,b). A similar 
calculation was carried out for the TR flow (Fig. 26a), where 
Cg = 2 also gives excellent agreement. Note that the predicted 
length scale variations model the integral scale changes as 
measured by CHC (Fig 25b). It therefore appears that a satis- 
factory model equation for the dissipation history in homogeneous 
flows is 

ft . . 4 4 + 2 4 (9.6) 

q q 

Next I considered the R^. equations with the objective of 
obtaining a model that, with (9.11) and (9.6), correctly predicts 
the measured R. . . Closure assumptions are required for the 
pressure strain and dissipation term. In all calculations I took 
(6.5) and (8.11) in (6.1), and hence wrote 


dR. 


dt 


i au. au i 

t 1 - ~ R ik 33^ - R jk + p ij 


3 e6 U 


(9.7a) 


where the pressure-strain term is 

du, 8u . 

P ij = p <3x7 + ^ 

J -L 

I first took 


(9.7b) 


P u ■ c 3 € 6 ij - h/9 2 ) 


(9.8) 


The first test was for the strain-free portion of the TR flow, 
where the structure is relaxing towards isotropy. Calculations 
showed that C^ = 5 gives a good representation of the structure, 
energy, and production in this flow (Fig. 26b). 
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I then proceded to try (9-8) with = 5 in the straining 
regions of the TR and CHC flows, but was not satisfied with the 
energy predictions (see Figs. 25c, d, 26c, d) . It appears that 
some alteration in either the P^. or dissipation terms is 
required, and I chose to experiment with modifications in the 
P ij * A ground rule was that proposed modification could not 
alter what has already been systematically established. The 
forms investigated were 

R 

p ij - ( 5e + C 4<P - • J 4) (9.9a) 

Q. 

R 

P ij = 5€ (i 5 ij “ ~2^ + °4 q2s ij (9.9b) 

Si 

hj ■ * . Jij 

+ + B jAl + § 9 »ij) = V (9.9c) 

dU. dU. dU. dU. 

[R ikfe^ “ 33T> + R jkfer " 377 

■K J K X 

(9.9d) 

Note that P^ =0 in each case. Eq. (9.9a) is DH's form with 
slightly different constants. Equations (9*9b)-(9.9d) are 
suggested by the notion that interactions between the mean strain 
rate and fluctuation fields contribute to the pressure fluctua- 
tions. The closures (9.9a) and (9»9b) were unsatisfactory. For 
the TR flow (9.9c) with = 1/2 works very well (Fig. 26c, d), 

but it is not adequate for the CHC flow (Fig. 25c, d). Equation 
(9»9d) reduces to (9*9c) for irrotational mean flow, i.e. the TR 
flow, and with C^ = 1/2 and C^ = 1/4 (9*9d) predicts the CHC 

flow reasonably well (Fig. 25c, d). 


P ij 


- V + °5 
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It does seem clear that (9.7) is not adequate in flows with 
strain or shear. With the constants indicated, (9.9 d) is 


R. 


P iJ - 5* (-J 6 ±J - 4 1 ) + + R Jk S kl + 

, £>u, au„ au, au v 

+ ¥ tR ik ( 53^ ' 3xj) + R jk^ • 33^ (9-lo) 


which is probably better. Further development is needed, and 
(9.10) is offered here as an interim model. However, it is not 
clear that (9.10) is a model of P^j ; it could just as well be a 
model for its complement (see 8.18) as used by Lumley (1970) in 
( 8 . 16 ).' 
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10. Outlook for the Future 


It should not be long before simple boundary layer flows 
are routinely handled in industry by MVF prediction methods. 

These methods are easy to use, require a minimum of input data, 
and give results which are usually adequate for engineering 
purposes. MTE methods will become increasingly important to both 
engineers and scientists, for they afford the possibility of 
including at least some important effects missed by MFV methods. 
The debate over the gradient-diffusion vs. large-eddy-transport 
closures will continue, and both methods will probably continue 
to be used with nearly equal success. MRS methods will be 
explored from the scientific side, and porbably will not be used 
to any substantial degree in engineering work for some time to 
come. 


Considerable effort is likely to be expended on the develop- 
ment of length scale (or equivalent) equations, such as the 
dissipation equation discussed in Section 8. In this connection 
the two -point correlation 

= Ui(x)Uj(x +0 (10.1) 

could be used to advantage, either along the lines of Gawain and 
Pritchard (1970), or perhaps through a closure of its own 
dynamical equation (Hinze 1959)* 

This equation will of course involve six independent space 
variables, but by integration over the separations £ these 
variables could be removed. Then, one might assume the form of 

R i j ' say 

ft. . = R. . e~^At (10.2) 
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and carry out the integrations* thereby obtaining three addi- 
tional differential equations of the transport type relating 
the integral scales and the one-point correlations 

(turbulent stresses) R. . . Experimental calculations along 
these lines would be most interesting. 

The heavy computation approach (Deardorff 1970) might be 
used to numerically test the closure assumptions used in the 
simpler MRS and MTE models. Hopefully such computations will be 
documented in the future with this use in mind. 

The MRS closures will attract most interest for use where 
MTE methods fail. For example* in flows with rotation the 
coriolis terms enter the R. . equations, but drop out in the 

p J 

equation for R. . = q . Therefore* it probably will be 
essential to go to a MRS method to include rotation effects* which 
are of considerable importance in many practical engineering and 
geophysical problems. Other effects which have not yet been 
adequately modeled and for which MRS methods may offer some hope 
include additive drag reduction* ultra-high Reynolds numbers* 
separation* roughness* lateral and transverse curvature* and 
strong thermal processes that affect the hydrodynamic motions. 

We might also see the complex closure models used as the 
basis for computationally simpler integral methods. The success 
of integral ' methods of this type at TBLPC (MC*HR; TBLPC-1) should 
not be forgotten in the rush to use the full partial differential 
equations . 

A disappointing aspect of the current status is that very 
little use has been made of the substantial advances made over the 
past decade in our understanding of the structure of turbulent 
shear flows. We know that large eddy structures dominate such 
flows; only MTES methods recognize this at all, and then not 
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quantitatively; MTEN methods ignore it all together. MRS cal- 
culations have not used the large eddy idea at all. Indeed, the 
concept of a "transport theory" for turbulent correlations would 
seem antithetical to a large eddy view. We know that the wall 
region is dominated by a particular correlatable structure, yet 
none of the current prediction methods make any use of this fact. 
The structure of the outer region of boundary layers has been 
extensively studied recently, and entrainment of non-turbulent 
fluid through the turbulent interface (superlayer) is known to 
be a critical process in turbulent shear flows; no real utiliza- 
tion of this fact has been incorporated into any of the closure 
models. We know that the outer layer flow has a dual structure, 
intermittently consisting of turbulent and non-turbulent regions 
with considerably different character. Yet all calculation 
methods are based on averages taken over long periods of time, 
averages that wash out this essential feature of the flow. Some 
believe a similar duality exists in the wall region; might not 
this also be incorporated? 

In short, it seems that too much attention has been paid to 
the numerical aspects of the computations; indeed, the difficulty 
of a first encounter with complex differencing schemes has made 
this necessary. But now we should begin a concerted effort 
to bring the new physical information into the turbulent flow 
computation methods, and we look for a better situation ten years 
hence . 

Finally, for those who cannot or choose not to take up the 
computation game, perhaps some fresh thinking at the fundamental 
level might be fruitful. For example, we all want to compute 
turbulent flows. Just what is it that we are working so hard to 
compute? What is the operational definition of turbulence? 
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